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1. Introduction 

Characteristic formulations of the Einstein equations have proven to be an important 
tool for numerical relativity. Most recently, they have been employed for the practical 
problem of measuring gravitational waves in a gauge-invariant and unambiguous way 
from numerically evolved spacetimes of binary black hole mergers [2 [21 IS1 E] , rotating 
stellar core collapse [5], and collapsar formation [6]. The technique, called Cauchy- 
characteristic extraction (CCE) (see [7] for a review) takes metric boundary data 
produced by a 3+1 evolution on a worldtube T of finite radius, and uses null-cone 
evolutions of the Einstein equations to transport that metric data to future null 
infinity, JJ + , the conformal outer boundary of spacetime where gravitational radiation 
is invariantly defined and interpreted in the Bondi gauge (HEH!!]]. Thus, this technique 
removes the influence of near-zone and coordinate effects [lj [3J [4), [5] . 

The null formulation is extremely efficient for evolving fields in the wave zone, 
where the null coordinates are well-behaved and caustics along geodesies are unlikely 
to be an issue. On relatively small computational grids, by comparison with standard 
3+1 methods, it is possible to achieve an accuracy which has proven to be sufficient 
for practical applications (gravitational wave measurement from compact bodies). In 
the case of characteristic gravitational-wave extraction, the inner boundary data for 
characteristic evolution is constructed on a worldtube at some distance form the source 
where the curvature gradients are already rather small compared to those close to 
the source. Thus, characteristic extraction requires comparatively little numerical 
resolution and is therefore less computationally demanding than 3+1 evolutions in 
the near-zone of a dynamical source. 

Although the computational effort is smaller, it is still non-trivial. To yield 
sufficient accuracy, for example to extract the gravitational radiation from a binary 
black hole evolution, the characteristic computation by current methods requires 
several days on up to a dozen processors on a workstation or small cluster to 
complete a 3000M time-series encompassing a dozen orbits of the binary (where 
M is the dimensionless mass of the spacetime). The application of higher-order 
discretization schemes, on the other hand, may deliver sufficient accuracy at much 
lower computational cost so that the additional effort of characteristic extraction 
would become negligible. If the characteristic code could be run concurrently with the 
underlying 3+1 simulation, this would allow for on-the-fly extraction of waveforms, 
as well as the interesting potential of using characteristic methods to provide exact 
non-linear boundary conditions for 3+1 codes. 

Characteristic codes have a long history in numerical relativity. A prominent 
result was the first stable dynamical evolution of a black hole spacetime in three 
spatial dimensions achieved by the Pittsburgh group [TT] [T2]. Since then, the 
Pittsburgh null code (or PITTNullCode) has become the main building block for 
current implementations of characteristic extraction used in numerical relativity 
simulations [IJ [H [31 [H [5] , and is now part of the publicly available Einstein Toolkit 
[13] . The code employs a single- null coordinate system, and is formulated in terms of 
spin- weighted variables that are related to the original variables defined by Bondi and 
collaborators 8, 9, 10 . It is built around a numerical scheme based on points located 
at the corners of null parallelograms, which was originally shown to be stable [H] in 
the context of the scalar wave equations. The PITTNullCode in its original form is 
2nd-order accurate in space and time discretization, meaning that the error decreases 
as 0(1/N 2 ) as the resolution N is increased. 
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Over the years, there have been several improvements to the original algorithm 
as implemented in [TTJ [T2] . Stereographic coordinates were replaced by more uniform 
angular grids |15l 116) . Fourth-order accurate angular derivatives have also been 
introduced, though retaining the 2nd-order parallelogram-based radial integration 
scheme [T51 15] . 

An alternative to the 2nd-order accurate null parallelogram was attempted by 
Bartnik and Norton [T71 [18]. They developed an algorithm based on a null quasi- 
spherical gauge, using method-of-lines integrators and 4th-order in time accuracy. 
They also used a spherical harmonic decomposition of variables on angular shells, and 
thus their code was pseudo-spectral in the angular directions. Their formulation of the 
null evolution equations (and coordinates) lead to certain numerical complications, 
including the need for high-accuracy interpolation operators to compute radial 
derivatives, and an elliptic gauge equation. Ultimately, the code did not demonstrate 
long-term stability. 

In this paper, we present a new high-order integration scheme for characteristic 
evolutions in general relativity. The evolution equations are written in the Bondi 
form, following the prescription of 12J. Time and outward radial integrations are 
performed by method-of-lines schemes. In particular, for the time integration, we use a 
classical 4th-order Runge-Kutta integrator. For the radial direction, we use a modified 
Adams-Moulton multi-step method. A multistep method is required by the lack of 
information at points between grid spacings, which is needed by the intermediate steps 
of the Runge-Kutta schemes. A similar method, for the case of 2nd-order accuracy 
and axisymmetry, was previously used in [19]. In the present context of full three- 
dimensional characteristic evolutions, we additionally need to discretize the angular 
direction. For this, we use spectral expansions in terms of spin-weighted real-valued 
spherical harmonics. The equations are solved on an angular collocation grid using the 
pseudo-spectral method. Radial and time integrations are performed on the angular 
spectral coefficients of the evolution variables. 

By referring to a simplified linear model problem in Section [2] we argue 
analytically that the proposed method is stable. We review the Bondi evolution 
equations in Section [3] In Section [4] we describe the numerical methods, including 
time and radial integrators, and pseudo-spectral derivatives. Finally, we test a newly 
implemented three-dimensional high-order code, demonstrating the expected order of 
convergence and accuracy of the method. We find that the new scheme is significantly 
more efficient than the old scheme, reaching the same level of accuracy using only very 
few radial and angular points. 

2. Integration schemes and stability 

Characteristic evolution schemes based on null or double-null coordinates have a 
significantly different character than conventional 3+1 (Cauchy) evolutions. In 
3+1 evolutions, spacetime is foliated along a timelike vector field t a by spacelike 
hypersurfaces S. In so-called 2+1+1 characteristic evolutions (the characteristic initial 
boundary value problem [20]), spacetime is foliated along a timelike vector field by 
null hypersurfaces, which are characteristic surfaces of the Einstein field equations. In 
this section, we describe a simplified model characteristic problem which exhibits the 
important features of the Einstein system, which we outline explicitly in Section [3] 
The characteristic method solves for the values of a field, J, which obeys a 
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hyperbolic equation of the form 

Here u = t — r is a retarded time coordinate labeling individual null slices, r is a radial 
coordinate along each null slice, and the index A = (9,<fi) labels angular coordinates. 
We allow the function F to depend on </, and its partial derivatives (to 2nd-order in 
r), which we label by subscripts 

d 2 J 

dudr ^ ' ur - [ ) 

We consider the problem on a domain bounded on the interior by a timelike worldtubc, 
r, at some finite areal radius R = R(&, <p) from the centre of our coordinate system, and 
on the exterior by S + ■ Boundary data consists of the variables required to evaluate 
F on r, as well as initial data for J along a single null cone at u = uq. 

We introduce an intermediate variable $ = J u , which allows us to recast Eq. ([IJ 
as a pair of lst-order equations 

$ r = F(J, J tr , J irr , J iA ) , (3a) 

Ju = $ • (36) 



Equation (3a) does not involve time derivatives J u and can thus be solved on a 
u = constant null slice by radial integration. Given data for J on a u = constant 
slice either from the last time step or by appropriate initial data, we propagate the 
system forward in time by first solving Eq. (3a) along the radial direction, and then 
integrating Eq. (3&) forward in time to determine J on the next slice. 

The schemes which we use for both radial and null integrations fall within the 
broad class of method-of-lines integrators for partial differential equations. These 
schemes assume that we have been able to evaluate the RHS at a point, so that 
standard ordinary differential equation (ODE) solvers can be applied to evolve the 
functions either forward in time, u, or (in the case of $) outward in r. 

For the time direction, we can use an explicit integrator, such as a standard 4th- 
order Runge-Kutta scheme. Recall, however, that in taking the solution from u to 
u + Am, the Runge-Kutta method involves calculating the right-hand side (RHS) at 
a number of intermediate steps. In the case of the J integration the RHS is $, and 
thus we need to compute $ via radial integration of Eq. (3a) at intermediate substeps 
between our timesteps of size Au. 

Applying such a scheme in the radial direction for $ is more problematic. In 
that case, the RHS is given by F, which is a function of data which is only known 
on discrete spheres (separated by a fixed distance, Ar), and so we cannot evaluate 
the intermediate substeps required by a Runge-Kutta-type integrator. Alternatively, 
we could make use of a multi-step algorithm, such as Adams-Bashforth or Adams- 
Moulton. These methods evaluate the RHS over some number of previous points, 
which in the case of the radial integration correspond to a set of the radial spheres on 
which J and its derivatives can been evaluated. 

We have examined a number of numerical schemes for carrying out the radial 



integration required by Eq. (3a) for the Einstein system, but almost universally 
found them to be unstable in empirical tests. To investigate the stability of different 
numerical methods, we turn to a simpler model which embodies the main features of 
the Einstein system, on which we can carry out a von Neumann analysis. 

By setting individual terms in the Einstein equations to zero and examining 
the subsequent numerical evolution, we came to the conclusion that the key terms 
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determining the stability are those which involve the variables J r and J rr . Thus 
we constructed two simplified systems which consist exclusively of these terms in the 
radial integration of a variable $ = J u . That is, we considered the systems 

$ r = Jr , J,u = $, (4) 

and 

$ r = J,rr , J,u = ®, (5) 

individually. In practice, we replace the areal radius r by a compactificd coordinate 
x, defined by 

x = — *— (6) 
rr + r 

where rr is a non-zero parameter corresponding to the radius of the inner boundary 



(see Eq. ( 25 ) , below) , and have found this transformation to be important in the 
stability analysis. In terms of x, Eqs. Q and |5]) become 

<E,x+ (r + rr)2 «D= rr J x , J U = <J>, (7) 
rrr r(r + rr) 

and 

(r -\- rrl rr 

*,* + - —^ = ^r^^2 J ^' J,u = ®, (8) 

rrr 2(r + rr) 

respectively. 

The von Neumann analysis corresponds to assuming the following form for the 
model variables: 

J = e wu e ikx , ^ = EJ, {9a) 

kh 

Ax = h, Au = fxh, v= ~2' *- 9 ^ 

The J u equations are evolved using a 4th-order Runge-Kutta integration, for which 
the general stability analysis is quite involved. However, to leading order in h, it is 
the same as the Euler method (and this was also found for other evolution algorithms 
such as the Adams-Bashforth methods). Thus, 

e w » h = 1 + [ihE, (10) 

and the stability is determined by the sign of 3t(E); 

> 0, Unstable, 

$t(E) { = : Stability unknown, (11) 

< : Stable. 

The quantity E is determined by the particular finite difference algorithm used to 
evaluate Eq. ([7| or Eq. pi), followed by substitution of Eqs. (9a 96). We investigated 



both second and 4th-order explicit (Adams-Bashforth) and implicit (Adams-Moulton) 
multi-step methods, with finite differences evaluated using centred, forward and 
backward methods of the appropriate accuracy. The calculations are somewhat 
lengthy and were done using a computer algebra script (see the file vN_comp3 .map 
in the online supplement). 

We found that the case Eq. ^ involving J jX always leads to values of v for which 
di(E) > so that the system is always unstable, independent of the integration scheme 
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for However, there are cases for which Eq. ([8]), involving J tXX , is stable. Since a 
second derivative term in the stability analysis is divided by h 2 compared to division 
by h for the first derivative term, to leading order in h the stability of the second 
derivative term dominates that of the first derivative term. 

The stability analysis of Eq. ^ indicates that among the tested methods, only 
the following cases are stable: 

• 2nd-order, Adams-Moulton, centred differences; 

• 2nd-order, Adams-Bashforth, forward differences; 

• 2nd-order, Adams-Moulton, forward differences; 

• 4th-order Adams-Moulton, forward differences; 

where by "forward" difference operator we indicate that radial derivatives of J are 
calculated using a stencil which involves points in the positive radial direction. We 
were able to confirm these results by empirical tests with the simplified system, Eq. 
Thus, for the full Einstein equations, described in the next section, we implemented a 
scheme in which radial integrations are carried out using a 4th-order Adams-Moulton 
method with upwinded radial derivatives. 

3. The Einstein equations in the Bondi-Sachs framework 

3.1. Coordinates 

The Bondi formulation writes the Einstein equations in terms of a null foliation of 
an asymptotically flat spacetime. We introduce coordinates y a — (r,y A ,u). The 
coordinate r is a radial surface area coordinate, and u = t — r is a retarded time 
coordinate which replaces the time, t, of 3+1 formulations. The y A are angular 
coordinates, labeled by uppercase indices that take the values 1 and 2. The coordinates 
(y A ,u) label individual null geodesies extending from a worldtube r = S 2 x K. The 
worldtube is chosen so that r = constant on T. In these coordinates, the general 
spacetime line element is 

ds 2 = - (e 2 ^- - r 2 h AB U A U B ) du 2 - 2e 2 ^ dudr 

-2r 2 h AB U B dudy A + r 2 h AB dy A dy B , (12) 

where h AB satisfies 

h AB h BC = 5 A c, det(h AB ) = det(q AB ) , (13) 

and q AB is the unit sphere metric. 

The spacetime is described by V, fi, U A , and h AB , which are functions of the 
coordinates. It is convenient to write quantities in terms of spin-weighted scalars in 
order to remove explicit angular tensor components. This simplifies the expression of 
the field equations in a way that is independent of the choice of angular coordinates. 
To this end, we introduce a complex dyad q A satisfying q A = q AB q B , q A q B = 0, 
q A q A = 2. [§j By projecting the angular variables onto this dyad, we define the 
complex valued scalars 

J = X - h AB q A q B , K=^h AB q A q B , U = U A q A , (14) 

§ An explicit form for q A will not be needed here, since we will represent angular dependence in terms 
of spin- weighted spherical harmonic basis functions constructed using a particular dyad representation 
adapted to the choice of angular coordinates. 
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of spin weights 2, 0, and 1, respectively. The components of hAB are uniquely 



determined by J due to the determinant condition, Eq. (13), thus fixing K as & 
function of J via 

K = VT+ J J . (15) 

Corresponding to the complex dyad, we introduce complex angular covariant 
differential operators 3 and <3 which maintain the property of spin-weight when 
acting on a scalar $ of spin- weight s [21]. The action of the 8 and 8 operators is 
restricted to transformations of our spin-weighted spherical harmonic basis functions 



(see Section 4.3 ) 



Before writing out the Einstein equations, we note that it is convenient to 
introduce an additional intermediate variable defined by 

Q:=r 2 e- 2 ^h AB U B r q A . (16) 

This spin-weight 1 variable, which is the first radial derivatives of U, will allow us to 
write the equations in lst-order form. Also, we re-express V in terms of a new variable 

V - r 

W:=—, (17) 
which has a regular limit as r — > oo. 

3.2. Einstein equations 

In Bondi coordinates, the vacuum Einstein equations 

Rab = (18) 

give rise to a hierarchy of equations which we can characterize as (i) hypersurface 
equations, (ii) evolution equations, and (hi) constraints. 

Hypersurface equations do not depend on u-derivatives, and thus can be evaluated 
within a u = constant slice. They are determined by the components R rr , RrAQ A , 
and RABh AB and lead to the following hierarchy of equations: 

0,r = , (19a) 
(r 2 Q), r = -r 2 (BJ + dK), r + 2r 4 B(r- 2 (3) r +N Q , (196) 

U r = r~ 2 e 2fi Q + N v , (19c) 
{r 2 W), r = ^TZ-l-e^me^ 

+ ^r- 2 (r 4 (W + W)) r + N w , (19d) 
where the Ricci scalar is given explicitly by 

n = 2K- + 1{B 2 j + g 2 J) + -r^(§J3J - B,mJ) , (20) 

2 4a 



and Np , Nq , Njj , and Ny/ are non- linear aspherical terms given explicitly in Appendix 
|A"} The equations are solved in succession, assuming available data J on a u — constant 
slice and constraint satisfying inner boundary data at the worldtube T for each of the 
hypersurface variables. This allows us to solve for /3, which in turn provides data for 
the equation for Q. Given j3 and Q, we can then solve for U, and finally for W. 
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The RabQ Q component of the Einstein equations determines the evolution 
equation for J: 

2 (^),„r-(( 1 + ^) ( r J),r) r = 

- r^ 1 (r 2 W) r + 2r- 1 e /3 g 2 e /3 - (rW) , J + Nj , (21) 
where the non-linear aspherical terms have been gathered in the quantity Nj (specified 



Appendix A I . We introduce an intermediate variable 
$ := J,u ■ 



(22) 



In terms of $, Eq. (21 ) becomes a new hypersurface equation 

(fi. + rW) (rJ), r ) r = 
~ l (r 2 W) r + 2r~ 1 e' 3 o i V - (rW) r J 



2 (r$) 

A ./ • (23) 

which is integrated radially from T using known values of the hypersurface variables 



determined in Eq. ( 19a|-( 19d). Then, J is determined by a timelike integration of 

J u = $- (24) 

We have expressed the Einstein system in a form analogous to the simplified model 
described in Eq. The source for 4> is complicated, but determined entirely by 
radial integration. Note the presence of the J rr in the second term of Eq. (23 1, which 



we have highlighted in Section [2] as key to determining the stability of numerical 
evolution schemes. 

Finally, we take advantage of the nature of null geodesies in asymptotically fiat 
spacetimes to compactify the radial direction so that is a boundary point of a 
closed domain. We replace the areal radius r by a new coordinate x via the invertible 
coordinate transformation 

x(r) = — T —, r{x) = r r ^—, (25) 

where rr is a constant, which we choose to be the radius of the worldtube, T. In this 
coordinate, the equations have a regular limit as x — > 1, and furthermore, we are able 
to set terms of order l/r n for n = 1,2, . . ., to zero at J + (see Section 3.3). Throughout 
the domain, derivatives are evaluated numerically in terms of the new coordinate, x, 
and then transformed into r-derivatives using the standard Jacobian transformations. 



For Eq. (25 1, these are 
dx 

dr (r 



r r ) 2 



d x 
dr 2 



2r r 
(r + r r ) 3 



(26) 



3.3. Form of the equations at J + 

In the limit of r — > oo, corresponding to i = 1, the numerical treatment of the 
equations requires special care. The problematic terms are those involving the 
Jacobian and the coordinate function r(x), which are not regular as r — s- oo. If 
the coordinate function and the Jacobian are explicitly inserted into the equations, 
using their form given by Eq. ( 25 1 and Eq. (p6| , respectively, divergent terms are seen 



to cancel. However, since we do not explicitly impose a specific compactification — 
we have formulated the problem in terms of generic Jacobians rather than the specific 
formulas of Eqs. ( 25 ) and ( 26 ) — we need to be careful to avoid irregular terms at J + . 
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It is sufficient to require that in the limit 



oo, the compactificd coordinate 



transformation and its Jacobian approach the explicit forms given in Eq. ( 26 ) and 

(27) 



Eq. (25), respectively. For instance, consider the equation 
(r 2 {x) d ° 



U 



dr 



)">(< 



r e 



According to Eq. (26) and Eq. (25), we have 



(r 2 (x) 



dx\-! 
~dr) 



as r 



oo. Hence, at J 
U 



r- x e^ 



(q 



+ rV 



2 ?N r 



r^e 2 ? 



(kq 



JQ 



(28) 



(29) 



We proceed in a similar manner for the other hypersurface equations. The specific 
form of the equations at is given in Appendix B Note additionally that for Q, 
W and $ it is possible to directly evaluate the respective quantity without radial 
integration at . For instance, as r — > oo, 



Q = -2573. 



(30) 



4. Numerical methods 

Discrete representation of the evolution variables 

The evolution algorithm is a hybrid of finite-difference (for radial and time integration) 
and pseudo-spectral (for angular directions) methods. In the compactified radial 
direction x, fields $ are evaluated on a uniform grid of N x points, i = 0, . . . , N x — 1 
at points x € [a;; n , . . . , 1]. The inner coordinate radius, x in = Xi=q, is that of the world- 
tube r, where we need to specify appropriate boundary data at any given time u to 
carry out a radially outward hypersurface integration. As we will see in Section |4.4| 
our radial integration scheme actually requires 3 radial points to start the algorithm. 
We therefore need to provide boundary data on the first i — 0, 1,2 radial points so 
that our worldtube T spans three radial points. Boundary data is required for 



Q, U,W,&\ for all 



i/i=0,l,2 



eri. 



(31) 



The outer boundary of the compactified radial grid is placed at the outermost gridpoint 
i = N x — 1 corresponding to future null infinity J + . 

At each radial point Xi, we represent angular dependence as a spectral expansion 
in terms of real- valued spin-weighted spherical harmonics, according to 



00 ra—-\-t 

EE 

ll=s m=-l 



®£m(Xi) s Zl m {y A ) 



(32) 



where the s Ze m are spin s real- valued spherical harmonics, defined in terms of the 
standard s Ye m basis [22] by 

1 

Zt m = { a Yi m , m=0 , (33) 



{ s Ytm + (~1)™ sYt-m) 



m > , 



V2 



(( — iy n s Y( m — s Yt-m) 



m < . 
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We use the s Zi m to accommodate the real- valued spin-0 quantities, which naturally 
yield real valued coefficients. 

We store a finite number of harmonic coefficients for each variable, terminating 
the sum according to the maximum number of measurable gravitational wave modes 
contained in the solution. For our test case with linearized solutions as discussed in 
Section [5] this is € max = 3. For the realistic case of current binary black hole merger 
simulations the number of resolved gravitational-wave modes in the 3+1 evolution is 
typically £ max ~ 8, beyond which their amplitude is below the level of discretization 
error. Although this case is not considered here and is left for future work, we do 
present results of a stability test in which £ max = 8. We store the spectral coefficients 
of the expansion of each evolution variable at each point of the radial grid. 

Radial and time integration is performed entirely on the spherical harmonic 



coefficients of the evolution variables, with the one exception being Eq. (56), to be 
introduced below. This equation contains non-linear terms of the form a $ where a 
and <£> are both functions of angular coordinates y A . To compute non-linear terms 



occurring either in Eq. ( 56 1 or in the RHS for a given hypersurface equation, we 
first need to recompose the involved variables on a collocation grid. After the terms 
have been evaluated on the collocation grid, we decompose them back into real-valued 
spin- weighted spherical harmonics. 

We construct a pseudo-spectral collocation grid for spherical harmonics by 
defining a set of grid points on a spherical shell S 2 using a y A — (9, <p) spherical- 
polar coordinate system with constant grid spacing in 9 and cf> direction 

9j,<f>k) = ( 7r ^ t ' 27r j^) : ^ fc e^0<J<^,0<fc<A^J . (34) 

Recomposing quantities on the collocation grid is easily done by evaluating the sum 
of the spherical harmonic expansion via 

f(Qj,<j>k) = / , fimsZimjOj, <f>k) ■ (35) 

tm 

Decomposing a quantity into spherical harmonics requires surface integration over 
S 2 . The expansion coefficients are computed according to a discrete version of 



hm = / dnf(9,q>) s Z em (6,cj)), (36) 
Jn 

> d8 d(j> is the surface element on the collocation grid. A numerical 
integration algorithm which is exact for spherical harmonics up to order (I, m) is given 
by Gauss-Chebyshev quadratures using (Ng, N<p) = (2(1 + 1), 2(0. + 1)) points on S 2 ( 
e.g. [23]). This algorithm makes use of coordinate dependent weights Wj in 9 direction. 
In <j), the weights are simply 1 since cj> is a periodic coordinate direction. Since we use 
equally spaced points in 9 (equivalent to Chebyshev nodes in x = cos 9), the weights 
in 9 direction are [53] 

4 ^A" 1 1 



e=o 



E ^jsin^ + l)^). (37) 



The integral Eq. ( 36 1 reduces to 

N e N, 



^ m = F¥"I!E f(0j,</)k)sZem(6j,^k) sinOj w 3 , (38) 

r j=0 fc=0 
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provided we have 

(JV fl ,JV*)<(2(* + l),2(*+l)) (39) 

points on S 2 . 

We have thus established an exact mapping between the representation in terms 
of spherical harmonic coefficients and the representation on the collocation grid. To 
speed up the computation, we precompute the s Zi m , and the product Wj sin^-. 

As an example for our algorithm, consider the linearized version of the 
hypersurface equation for U, 

U x =(r 2 (x)^y\^Q. (40) 

We recompose f3 and Q on the sphere using their spectral expansion coefficients 
in order to carry out the required multiplications. The complete procedure for 



integrating Eq. (40) can be summarized as follows. 

1. Loop over radial points. On each radial shell xi\ 

i. Define grid points on a spherical shell S 2 according to Eq. (34). 



ii. Recompose 

Pfa^fa) =J2P £m (xi) oZimiOj^k), (41) 
tm 

and 

Q(x l ,6 J7 cj )k )=J2Q em (x l ) iZUOj^k). (42) 

t III 

iii. Loop over all angular points 6j and (f>k- For each angular point, compute 



U x (e 3 ,<j) k ) using Eq. (|40j. 
iv. Decompose U x (0j,(f>k) via 

7T 2tt Ne N * 

{Utm),x{ X i) = ^^y^y^ U,x(0;j, 0k) S Ze m (0j,(j>k) SmO.jW.j . (43) 
" V j=0 fc=0 

2. Radially integrate (Ui m ) iX (xi) to obtain L^ m (x,). 
The last step, radial integration, is described in more detail in Section |4~4) 

4-2. Radial derivatives and dissipation 

Radial derivatives of all hypersurface quantities are generally obtained from the RHS of 
their corresponding radial ODE integrations and hence do not need to be recomputed 
by means of finite difference operators. However, the metric variable J (and also 
K) itself is not directly obtained via radial integration and hence must be computed 
everywhere. We approximate J x and J tXX by means of finite difference operators of 



4th-order. The radial derivative of K can be obtained by using Eq. (15). 

According to stability analysis and empirical findings (see Section [2j, we apply 

fully side- winded derivatives with the stencil points in the direction of J + . We use 
4th-order first and second derivatives 

dfi = ^ + 4/ i+1 - 3/ i+2 + - , (44) 

«2* 1 f 15 , 77 . 107 , 61 . 5 . \ 

~Kx 2 \ ~4~ 6~ ~Q~J i + 2 _ Wi+3 + ~^Ji+i _ I > l 4o 7 
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where Air is the grid spacing in the compactificd radial coordinate direction. 
Close to J + when i > N x — 5, we switch to 4th-order centred stencils 

dh = ^ (+^/-2 ~ \h-2 + pi+i - ^/ l+2 ) , (46) 

= Aa?2 (~ V2^~ 2 + 3^ i_1 ~ 2^ + 3^ +1 _ 12 ^ i+2 ) ' 

and when i > N x — 3, we switch to side-winded stencils pointing towards the inner 
boundary 

dfi = ^ O^fi - ifi-i + - |/<-3 + , (48) 

«2, 1 ^!5, 77. 107 , 61 . 5 . \ .... 

9 /* = ^2 " T^ 1 + "e"^ -2 " 13/i " 3 + 12 /l ~ 4 " 6 fi ~ B ) ' ( } 

In addition, we apply a numerical dissipation operator to J. We use a 5th-order 
Kreiss-Oliger dissipation operator of the form 

D h = tSx" (/*-3 - 6 ^-2 + 15 /*-i - 20/, + 15/ i+ i - 6/ i+2 + f i+3 ) , (50) 

where edi SS controls the strength of the applied dissipation operator D. At the outer 
boundary (at where we do not have enough points to compute the dissipation 

operator, we use one-sided stencil derived for an overall 4th-order accurate summation- 
by-parts (SBP) operator (though we do not make explicit use of the SBP property). 
The particular stencil coefficients are derived in [25] . We explicitly state the stencil 
coefficients in |Appendix C| 

Empirical tests have shown that radial dissipation applied to J is crucial to 
improve the stability properties of our scheme. 



4-3. Angular derivatives 

Numerical derivatives in the angular direction are obtained via analytic angular 
derivatives of the spin-weighted real-valued spherical harmonic spectral basis 
functions. The action of the <3 derivative on the real-valued spherical harmonics is 
given by [52] 

3 sZtm = y/(e + S + l)(e-s) s+1 Z im . (51) 

The action of 9 on a spectrally expanded function is 

3/ = -7W(* + s + W - s ) s+iZi m ■ (52) 

Similarly, the action of 3 is given by 

8 ,Z tm = -y/W-s + !)(£ + s) s ^Z lm . (53) 



4-4- Radial integration 

The hypersurface equations are integrated in the radial direction using a multistep 
method. The classes of methods that we have studied for this problem are either 
the explicit Adams-Bashforth methods, the implicit Adams-Moulton methods, and 
a combination of the two in the form of a predictor-corrector scheme. However, as 
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discussed in Section [2j at 4th-order, explicit methods are unstable for our particular 
set of equations. 

Thus, the radial integration uses a fully implicit method. Fortunately, the 
Einstein equations in Bondi-Sachs form are particularly convenient for this purpose, 
as they form a hierarchy (as outlined in Section [3]) with only one unknown function in 
each equation. Furthermore, the equations are linear in this unknown. Schematically, 
we write each equation in the form 

dy 



dx 



yg(x) = f(x), 



(54) 



and the 4th-order fully implicit Adams-Moulton scheme can be written in explicit 
form as 



2/4+1 i 1 + \ h9i 



19 

24 



(fi - Vi9i) 



24 (/i-i 



1 



Vi-igi-i) + ^(/i-2 



Vi-29i-2)j- 



(55) 



We again note that the quantities we work with are the spherical harmonic coefficients 
of each variable, which are which are functions of the compactified radius Xi, according 
to the procedure outlined at the end of Section |4.1| 

The radial $ integration requires special treatment due to the nonlinear term, 
Nj in Eq. (21 ) (given explicitly in Appendix A). It is a function of both J M and J M , 
and thus, according to Eq. (22), both $ and <i>. We write Eq. (21) in the form 



where 



b = 



(r(x) 
J 



dx\ - 1 
dr) 



J 



2K 



2K 

{JK X - J X K). 



(JK <X - J X K) 



(56) 



(57) 



and i?$ contains the remaining terms but does not involve $ or $ or their derivatives. 
This is integrated using the scheme of Eq. ( 55 ) to obtain an equation that involves 
both and Qi+i. Taking the complex conjugate leads to a second equation in the 
two unknowns, and solving the system gives 



$4 



1 



Zh. 



= TAl 



3/i 

*" IT 



1 



3h 



Cti+l 

3h 



3h 



+ i 



a i+i ) — Ti—bi+i 



(58) 



where 



Ti = $ 4 



3h 



19h 

5h 
24 
h 



i - ®i-ib 



•i-i) 



24 



(59) 



The scheme above does not allow us to work directly with the spherical harmonic 
coefficients of $, K, J, K x , J x , and due to the non-linear terms <&a and l>6 
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which must be evaluated on the collocation grid. To compute Eq. ( 59 1 , we therefore 
recompose K, J, K x , J x , and on the collocation grid defined by Eq. ([34] ) to 
perform the required multiplications. Having evaluated $ according to Eq. ( |59[ ) on 
the collocation grid, we decompose $ to obtain its spectral coefficients in terms of spin 
s = 2 real- valued spherical harmonics for each radial point a;,. 

Note that the radial integration schemes Eq. ( 55 ) and Eq. ( 59 1 both require data 
on 3 radial points to start the algorithm. These must be provided as boundary data 
on the worldtube T. 



4-5. Time integration 

The evolution equation for J has the form 

J,u = $ ■ (60) 

This equation can be straightforwardly integrated via a 4th-order Runge-Kutta scheme 
using the spectral coefficients of $. In addition, we add numerical dissipation to 



Eq. (60). To be explicit, we solve 

(Jt m ),u = $t m +DJt m , V£,m,i, (61) 



where D is a dissipation operator defined in Section |4T2| Since it is necessary to solve 
the hypersurface equations to obtain the $f m , the hypersurface equations must be 
solved for each intermediate Runge-Kutta step. 



4-6. Summary of algorithm 

1. Assume data for J in the form of spectral coefficients J lm {xi) at (intermediate) 
timestep t n for each radial shell Xi . If t n is the first timestep, the Je m (xi) are 
given by initial data. 

2. Compute Kz m (xi), as well as radial derivatives J l ™{xi) and J l ™(xi) from J(xi) 
by means of Eq. ( 44 ) and Eq. ( 45 1 , respectively. Since K is related non-linearly to 
J, we need to recompose J(xj) from J em (xi) to evaluate K(xi) on the collocation 
grid. Afterwards, we decompose K(xt) to obtain K (xi). 

3. Provide inner worldtube boundary data for /3, U, W and $ at (intermediate) 
timestep t n in terms of spectral coefficients on the first 3 radial points. 

4. Integrate hypersurface equations in the order (i) j3, (ii) Q, (iii) U, (iv) W, and 
(v) $ by using the steps described in Section [1] and Section 4.4 

5. Evaluate next Runge-Kutta step for J jU = $ to obtain J at next (intermediate) 
step t n+ i as described in Section 4.5 



4-7. Remarks on the computational implementation 

We have implemented a new code within the Cactus computational toolkit [2"6l 127] , 
The underlying grid array structures are provided by Carpet [281 I29j . Memory 
handling of the collocation grid, and recomposition/decomposition in terms of spin- 
weighted spherical harmonics is provided by SphericalSlice |30| . 

Since the memory consumption of the implemented code is rather low, and since 
the computational efficiency is high, we do not currently decompose the domain to 
distribute the work load across multiple processing units. We do, however, make use 
of multi-threading via OpenMP to enable faster processing on shared memory units. 
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In particular, we use multi-threading in the following two kinds of loops: (i) when 
looping over spectral coefficients to perform, for instance, radial integration for each 
separate mode, and (ii) when looping over points on the collocation grid to evaluate 
non- linear terms in the equations (such as RHS evaluation). Depending on the number 
of spectral coefficients and the number of available cores within one shared memory 
unit, the observed scaling can be close to the optimum (though we remark that we 
have done a rather limited number of tests on a compute node with up to 12 shared 
memory cores). 

The implemented code is designed such that it can run concurrently with our 
3+1 evolution code Llama |311 132j . This is important for future application in on- 
the-fly Cauchy-characteristic extraction, where the metric data will be transported to 
Sf + during Cauchy evolution without the need of an additional post-processing step 
(which is currently the case for the algorithm presented in [HG3). Furthermore, this 
is necessary for a future implementation of Cauchy characteristic matching[3Sl H] m 
which the characteristic evolution is used to provide on-the-fly boundary data for a 
Cauchy evolution. 



5. Results 



5.1. Linearized solutions 

To test the convergence of our numerical scheme on a dynamical spacetime, we use 
solutions to the linearized Einstein equations in Bondi-Sachs form on a Minkowski 
background (Section 4.3 of [34]). We write 

J Iin = y/(£-me+l)(e + 2) 2 Z lm n(J e (r) e™), 

U lin = ^JW+T) x Z lm R(U e (r) e wu ), 

/3 lin = Z im ^((3,e wu ) 7 

W lin = Z em ^(W e (r)e wu ), (62) 

where Jt(r), Ui(r), fit, Wi{r) are in general complex, and taking the real part leads 
to cos(z^u) and sm(vu) terms. The quantities /3 and W are real; while J and U are 
complex due to the terms (3 2 Zi m and (5Zi m , representing different terms in the angular 
part of the metric. We require a solution that is well-behaved at future null infinity. 
We find [15] , in the case I = 2, 

& = A>, 

24/3 + ZivCx - iv 3 C 2 Ci C 2 



Mr) 
Mr) 
W 2 (r) 



36 4r 12r 3 ' 

-24ii//3 + 3^ 2 Ci - v 4 C 2 2/3 C x ivC 2 C 2 

36 " + ~V + 2^2 + 3r 3 + 4>4' 

24w/3 - 3j/ 2 Ci + v^C 2 ZivCi - 6/3 - iv z C 2 
6 + ~ 3r 

v*C 2 ivC 2 C 2 

- — + -r + ^ (63) 

with the (complex) constants /3o, C\ and C 2 freely specifiable; and in the case t = 3 

A = A>, 

n = 60/3 + ZivC x + v 4 C 2 Ch_ _ ivCh _ 
3 ^ " 180 lOr 6r 3 4r 4 ' 
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U 3 (r) 



W 3 (r) 



-60iz//3 + Zv 2 C x - w 5 C 2 
180 

2/3 Ci 2^ 2 C 2 5ivC 2 C 2 
_l — ™ _| 1 1 _| 1 _| f. 

r 2r 2 3r^ 4r 4 
60ii//3 - 3^ 2 Ci + ii/ 5 C 2 ifCi - 2/3 + v A C 2 



15 3r 
i2^ 3 C 2 4u> 2 C 2 hvC 2 iC 2 



(64) 



We establish convergence by testing the evolution quantities against linearized 



solutions listed in Section 5.1 The linearized solution provides initial data for J on 
a null cone, as well as boundary data for j3, Q, U, W, and $ at the worldtube T. 
During evolution, we compute the error e in all evolved quantities by comparing with 
the linearized solution. Since the code solves the general nonlinear case whereas the 
exact solution satisfies the linearized Einstein equations, we expect e(J) to converge 
towards zero at the order of accuracy of the numerical scheme only in a regime in 
which |e(J)| is much larger than any nonlinear contribution. 

We have performed a number of test cases using (£, m) = (2, 2) linearized 
solutions, (£, m) = (3, 3) linearized solutions, and a superposition of both. In the 
latter case, we compute J via 

J lin = ^ »( Mr)e wu ) , (65) 

where ki = \)£{t + 1)(£ + 2). The remaining superposed linearized solutions for 

all other quantities are constructed in the same way, using the appropriate spin weight 
for the sZ£ m and the appropriate ^"-dependent coefficients m, respectively (compare 



Eq. (62)). 

The linearized solutions depend on free parameters C±, C 2 , f5 and v which we 
have tested for a range of different values. Note that the amplitudes C\, C 2 , and /3o 
must be of linear order (< 10~ 5 ). 

In all cases considered, we find better than 4th-order convergence until the error 
roughly reaches the square of the amplitude of the linearized solution, beyond which 
convergence deteriorates, as expected, due to the emergence of nonlinear behavior. 

As a particular example, we show convergence of a superposed {£, m) — (2,2) + 
(3, 3) solution with parameters 

C\ = 3 x 10~ 6 , C* 2 = 1 x 10~ 6 , (66) 
ft=ix 10~ 6 , u=1.0. (67) 

Fig.[l]plots the L 2 -norm of the error ||e(J)|| on two resolutions rO and rl (see Table]!]) 
scaled for 4th-order convergence. We define the L 2 -norm in terms of the sum over all 
modes and radial points by 




\J\\ = J^(hm(xi)) 2 . (68) 

The appropriate convergence scaling can be determined from the convergence rate 
defined in terms of the grid spacing Ax, 

<«•> 
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where p — 4 is the expected order of convergence. By doubling the resolution, we 
expect the higher resolution error, rl, to be smaller by a factor of 16 given our 4th- 
order accurate algorithm, i.e., we should get 

As shown in Fig. [T] this is indeed the case. We measure better than 4th-order 
convergence (see further below for a discussion). Furthermore, the evolution is still 
stable after T = 20000 M (corresponding to ~ 6400 cycles of the solution). 

The grid settings and parameters for this test are given in Tabic [T] In all cases, 
we apply radial dissipation of amplitude e<ji S s = 0.2. The inner boundary is located at 



Rr = 15M (corresponding to rr — 15, Eq. (25 1). The inner compactified coordinate 
radius x- ln is chosen such that the nominal grid (i.e., the grid excluding the 3 inner 
boundary points) starts at Xi = ^ = 0.36. 

In Fig. [2j we plot the time L2-norm of the error e( J), defined in terms of the sum 
of all modes on all radial points over all time steps t n 

1/11 = J > ' {Um{xut n )f. (71) 




We consider radial resolutions N x = [13,17,21,25,29,33,37,41,45,49,65] with 
appropriately adapted time resolutions. As the resolution is increased, the error drops 
as expected. Note that even on the coarsest radial grids, the code achieves an accuracy 
with a relative error better than e( J) f=a 10 -4 (the amplitude of our solution is on the 



order of 10 -6 ; compare Eq. (66 1). This is a significant improvement in efficiency over 
the original 2nd-order scheme developed in [12] including its advancements {15 , 3 . We 
also note that when further increasing the radial and time resolution, the error falls 
below ~ 10 -12 and we observe an expected drop in convergence order due to nonlinear 
effects. By fitting a line through the sampled error norms (red dashed line in Fig. {21), 
we can measure the convergence order. In the present case, we find ~ 5. Note that 
we have excluded the coarsest and the finest three resolutions from the fit. 

We have also checked convergence for each individual hypersurface equation by 
systematically setting all other quantities to their linearized solutions, and have found 
that while the equations for f3, Q, $ and J converge consistently at 4th-order, the 
equation for U and W both achieve measured convergence orders of ~ 5 at the given 
resolutions (the coefficient of the 5th-order error is apparently larger than the 4th- 
order term). These contributions lead to an overall 5th-order convergence for J. 
With sufficiently high resolution (N x > 60 points) the 5th-order error term in the 
individual equations for U and W diminishes sufficiently so that we observe the 4th- 
order convergence expected of the algorithm. As we increase the resolution further 
to the point that the error approaches the square of the amplitude of the linearized 
solution, we start to see the influence of nonlinear terms and the comparison with the 
linearized exact solution no longer holds. 



5.2. Pseudo random noise 

A strong test for stability involves injecting (pseudo) random noise into the evolution, 
via the initial and boundary data [351 • By this method, any exponentially growing 
error modes, if present, will be stimulated at a much higher amplitude than would 
naturally occur due to truncation or round-off error. 
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Table 1. Numerical grid settings used. N x denotes the number of radial points on 
the nominal compactified grid, / ra ax is the number of angular spectral coefficients 
used, No and N$ denote the number of angular collocation points in 9 and <j> 
direction, respectively, and Au is the time resolution. Note that our choice of the 
number of angular spectral coefficients and collocation points results in an exact 
representation of angular derivatives and functions in the case of a solution with 

e < 3. 





N x 




Au 


rO 


17 


3 8 8 


0.1 


rl 


33 


3 8 8 


0.05 



(£,m) = (3,3) + (2,2) solution — r 

rl xl6 




t[M] 



Figure 1. L2 norm Eq. \6S\ of the error in the evolution variable J of a 
superposed (£, m) = (2, 2) + (3, 3) linearized solution on two resolutions rO and rl 
(Table Y\\ . Given that the amplitude of J itself is on the order of 10 -6 (compare 
Eq. |66f ), we observe that the relative error at the given resolutions is < 10 -4 . 
The error in the high resolution rl is scaled for 4th-order convergence. We observe 
better than 4th-order convergence (see text for a discussion). The evolution is still 
stable after T = 20000M. This corresponds to ~ 6400 cycles in the solution. The 
inset shows a close up of the error for the last 50M of evolution. 



In this test, we add noise to all spectral coefficients of J on the initial null 
hypersurface, and also to all spectral coefficients of all remaining quantities that are 
needed at the inner boundary, the worldtubc T, at each timestep. We add noise of 
amplitude A no j S0 = 10~ 2 to a (£, m) — (3, 3) linearized solution with parameters given 
by Eq. ( [66] ). Note that the amplitude of the noise is 10,000 times stronger than the 
amplitude of the linearized solution itself, and thus of nonlinear scale. 

In Fig.|3j we show the norm Eq. ( 68 ) of the error e( J) over a period of T = 2500M, 
where we have chosen resolution rl as listed in Table [l] To allow for higher frequency 
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Figure 2. L2-norm over time Eq. l |71| of the error in the evolution variable 
e( J) of a superposed (£, m) = (2, 2) + (3, 3) linearized solution on multiple radial 
resolutions. As the radial resolution is increased starting from N x = 13 up to 
N x = 65 points, the error decreases at better than 4th-order. The two highlighted 
(red) markers correspond to resolutions rO and rl (TableQ, respectively. The red 
dashed line corresponds to a fit excluding the first point and the last four points. 
We measure a convergence order of ~ 5. When the error reaches roughly the 
square of the amplitude of the linearized solution (~ 10 — 12 ) where quadratic (and 
higher order) terms become important, the convergence order starts to deteriorate. 
This is expected since the linearized solutions only satisfy the Einstein equations 
to linear order. 



angular modes, in the solution we use spectral coefficients up to £ — 8, and increase 
the number of collocation points to (Ng,N^,) = (18,18) accordingly. We inject noise 
not only into the (£, m) = (3, 3) solution mode, but also into all other modes up to 
t = 8, which would otherwise be zero. To control the stabilitjfjj] of the scheme, we 
set the amount of dissipation to ediss = 0.5. The test demonstrates that the non- 
linear coupling of angular modes does not lead to unstable behavior over the observed 
timescale. 

6. Discussion 

We have developed and implemented a high-order algorithm for numerically integrat- 
ing the full non-linear Einstein equations along characteristic null hypersurfaces in 
the Bondi-Sachs framework. The implemented code evolves the full Einstein equa- 
tions in the wave zone of a compact body, remaining stable, convergent, and achiev- 
ing high accuracy with relatively low computational cost compared to the previous 
2nd-order algorithm. Radial integration is performed in terms of a modified Adams- 
Moulton scheme that is 4th-order accurate. Radial derivatives are computed in terms 
of 4th-order finite differences. Angular derivatives are computed in terms of spec- 
tral expansions of real-valued spin-weighted spherical harmonics. It is in principle 
straightforward to design an algorithm along the lines used here at even higher order 
of convergence than 4th-order. However, the stability properties of the algorithm are 
clearly dependent on its order, and achieving stability for an algorithm of higher order 
may be problematic. 

|| We have found that lower values of e may trigger an instability at . This can be cured by either 
using larger e everywhere, or by just increasing e at J7+ . 
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Figure 3. L2 norm Eq. | |68| of the error in the evolution variable J of 
a (i,m) = (3,3) linearized solution with pseudo random noise of amplitude 
j4noisc = 10~ 2 injected into modes I < 8 at the worldtube. The error remains 
bounded after a few tens of thousand iterations. This indicates that the scheme 
is stable even in the presence of strong noise. The inset shows a close-up of the 
error over a timcscale of 50M. 



The implemented algorithm is a first step towards more efficient gravitational- 
wave extraction algorithms via Cauchy-characteristic extraction. It is also a first step 
towards Cauchy-characteristic matching in which the characteristic evolution is used 
to provide on-the-fly boundary data for a 3+1 evolution. In view of this, we have 
designed our code such that it can run concurrently with our 3+1 evolution code. 
The next step will be the implementation of an improved algorithm for worldtube 
boundary data transformation that couple Cauchy and characteristic evolutions, and 
that is of higher than 2nd-order. 
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Appendix A. The nonlinear terms in the Einstein equations 



The nonlinear terms Np,NQ,Nu,Ny/ and Nj in Eq. (19a) through Eq. (21) were 
first presented in [12]. We repeat them here, but with a mis-print in Eq. (A3) of [12] 
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e 2f3 

N V = — (KQ-Q- JQ) , 
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(A.l) 
(A.2) 



N Q = r 2 (1 - K){dK <r + 5J r ) + Q(JJ. r ) + B{JK, r ) - J <r dK 



(A.3) 



N w = e 213 [ (1 - K){mp + (5(38/3) + ^ [ J(§f3) 2 + J (30) 



-^(d0(pK-dJ)+dp(pK-dJ)j + ±[ JO- i Jd 2 i 



-(2KU, r U ir + JU 2 r + JU 2 r ). 



(A.4) 



where 



Nj = Nj 1 +Nj 2 +Nj 3 +N M +Nj 5 +N,je+Nj r +-(P 1 +P 2 +P 3 +P 4 )(A.5) 



N. n 



N J2 



N 



,74 



^ ^K(dJd/3 + 2dKd/3 - 8 J9/3) + J(8 J§/3 - 2dKd/3) - JQJQ(3^ , 
1 



3 J(rU, r + 2U) + 3J(rU, r + 2U)\, 



N J3 = (1 - ^)(rSZ7 ;r + 2W) - J(rdU r + 2917), 



= T ^e- 2 ' j (k 2 U 2 t + 2JKU r U r + J 2 U 2 )j . 



N J5 = --J r (W + BU), 

N J6 = r^(UQJ + UfiJ)(JJ r - JJ r ) 

+ (JK ir - KJ r )UdJ - U(dJ r - 2KdKJ r + 2JdKK tr ) 
- C/(9J r - KSJJ >r + JSJK tr )j , 

N J7 = r(J, r K - JK, r ) (u(SJ - BK) + U{M - 5 J) 

+ K(W - W) + (JBU - JW)^j , 
Pi = r 2 (^(J.rK - JK r ) + J -^{J, r K - JK r )j - 8 (r + r 2 w) {3, 
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P 2 =e 



P 3 
Pa 



+ j(jrp + (spy) + j(pp + (spy) + (g ,mp + a jsp) 



{rW. r + 2W) + {rW. r + 2W) 



e~ 2fl {2KU ir U ir + JU\ + JU'i). 



(A.6) 



Appendix B. Regularized equations at J + 



At J + , the equations simplify when inserting Eq. (25) and Eq. (26) and taking the 



limit r — > oo. In particular, terms containing powers of 1/r, l/r™ vanish. Below, 
we give the equations evaluated at J + . 



P,x 

Q 

u x 

w 



where 



N 



0, 

-23/3, 

r^e^iKQ-JQ), 
dU + W 

2 ' 
- W + N J>J+ , 

UdJ + UdJ 



J.J+ 



+ (1 - K)W- 



J(W-W) 



(B.l) 
(B.2) 
(B.3) 

(B.4) 
(B.5) 

(B.6) 
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Appendix C. Dissipation operator stencil 

Close to the outer boundary where we do not have enough points to apply the 



standard centred-stencil Kreiss-Oliger dissipation operator Eq. (50), we use side- 
winded dissipation operator stencils derived for SBP operator Z?4_2 of [55]. This 
particular dissipation operator is defined via coefficients ay and In the interior of 
the domain the operator reads 



AijUj 



^diss 



qo Ui 



7 
3=1 



H+j) 



(C.l) 



where the qi are the coefficients given in Eq. (50). Near the outer boundary (on points 
i = N x — 5, N x — 1), the operator can be written as 

7 
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where the coefficients taken from [25] reported below. We note that ediss > 

selects the amount of dissipation and is usually of order unity. 

/ 



(<Hf) 



48 
17 


144 
17 


144 
17 


17 











144 
59 


480 
59 


576 
59 


288 
59 


48 
59 








144 
4.3 


576 
43 


912 

43 


720 
43 


288 
43 


48 
43 





48 
49 


288 
49 


720 

49 


960 
49 


720 

49 


288 
49 


48 
49 



(C.3) 



References 



[1] C. Reisswig, N. T. Bishop, D. Pollney, and B. Szilagyi. Unambiguous determination of 
gravitational waveforms from binary black hole mergers. Phys. Rev. Lett., 103:221101, 2009. 

[2] C. Reisswig, N. T. Bishop, D. Pollney, and B. Szilagyi. Characteristic extraction in numerical 
relativity: binary black hole merger waveforms at null infinity. Class. Quant. Grav., 
27:075014, 2010. 

[3] M.C. Babiuc, B. Szilagyi, J. Winicour, and Y. Zlochower. A Characteristic Extraction Tool for 

Gravitational Waveforms. Phys. Rev., D84:044057, 2011. 
[4] M.C. Babiuc, J. Winicour, and Y. Zlochower. Binary Black Hole Waveform Extraction at Null 

Infinity. Class. Quant. Grav., 28:134006, 2011. 
[5] C. Reisswig, CD. Ott, U. Sperhake, and E. Schnetter. Gravitational Wave Extraction in 

Simulations of Rotating Stellar Core Collapse. Phys. Rev., D83:064008, 2011. 
[6] C. D. Ott, C. Reisswig, E. Schnetter, E. O'Connor, U. Sperhake, F. Loftier, P. Diener, 

E. Abdikamalov, I. Hawke, and A. Burrows. Dynamics and Gravitational Wave Signature of 

Collapsar Formation. Phys. Rev. Lett., 106:161103, April 2011. 
[7] Jeffrey Winicour. Characteristic evolution and matching. Living Rev. Relativ., 8:10, 2005. 

[Online article]. 

[8] H. Bondi, M. G. J. van der Burg, and A. W. K. Metzner. Gravitational waves in general 
relativity VII. Waves from axi-symmetric isolated systems. Proc. Roy. Soc. A, 269:21-52, 
1962. 

[9] R. K. Sachs. Gravitational waves in general relativity. Proc. Roy. Soc. A, 270:103-126, 1962. 
[10] R. Penrose. Asymptotic properties of fields and spacetimes. Phys. Rev. Lett., 10:66-68, 1963. 
[11] N. T. Bishop, R. Gomez, L. Lehner, and J. Winicour. Cauchy-characteristic extraction in 

numerical relativity. Phys. Rev. D, 54:6153-6165, 1996. 
[12] Nigel T. Bishop, Roberto Gomez, Luis Lehner, Manoj Maharaj, and Jeffrey Winicour. High- 
powered gravitational news. Phys. Rev. D, 56(10):6298— 6309, 15 November 1997. 
[13] F. Loffter, J. Faber, E. Bentivegna, T. Bode, P. Diener, R. Haas, I. Hinder, B. C. Mundim, 

C. D. Ott, E. Schnetter, G. Allen, M. Campanelli, and P. Laguna. The Einstein Toolkit: a 

community computational infrastructure for relativistic astrophysics. Classical and Quantum 

Gravity, 29(11) :115001, 2012. 
[14] Roberto Gomez, Jeffrey Winicour, and Richard Isaacson. Evolution of scalar fields from 

characteristic data. J. Comput. Phys., 98:11 - 25, 1992. 
[15] Christian Reisswig, Nigel T. Bishop, Chi Wai Lai, Jonathan Thornburg, and Bela Szilagyi. 

Characteristic evolutions in numerical relativity using six angular patches. Class. Quantum 

Grav., 24:S327-S339, 2007. 
[16] Roberto Gomez, Willians Barreto, and Simonetta Frittelli. A framework for large-scale 

relativistic simulations in the characteristic approach. Phys. Rev., D76:124029, 2007. 
[17] Robert Bartnik. Einstein equations in the null quasispherical gauge. Class. Quant. Grav., 

14:2185-2194, 1997. 

[18] Robert A. Bartnik and Andrew H. Norton. Einstein equations in the null quasi-spherical gauge 

III: numerical algorithms, gr-qc/9904045, 1999. 
[19] Nigel T. Bishop, C. Clarke, and R. d'Inverno. Numerical relativity on a transputer array. Class. 

Quantum Grav., 7(2):L23-L27, February 1990. 
[20] J. M. Stewart. Advanced general relativity. Cambridge University Press, Cambridge, 1990. 
[21] Roberto Gomez, Luis Lehner, Philippos Papadopoulos, and Jeffrey Winicour. The eth formalism 

in numerical relativity. Class. Quantum Grav., 14(4):977-990, 1997. 
[22] J. N. Goldberg, A. J. MacFarlane, Ezra T. Newman, F. Rohrlich, and E. C. G. Sudarshan. 

Spin-s spherical harmonics and 5. J. Math. Phys., 8(11):2155— 2161, 1967. 



General relativistic null-cone evolutions with a high-order scheme 24 

[23] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes in C++ 
: the art of scientific computing. New York, 3rd edition, 2002. 

[24] James R. Driscoll and Dennis M. Healy, Jr. Computing fourier transforms and convolutions on 
the 2-sphere. Adv. Appl. Math., 15(2):202-250, 1994. 

[25] Peter Diener, Ernst Nils Dorband, Erik Schncttcr, and Manuel Tiglio. Optimized high-order 
derivative and dissipation operators satisfying summation by parts, and applications in three- 
dimensional multi-block evolutions. J. Sci. Comput., 32:109-145, 2007. 

[26] T. Goodale, G. Allen, G. Lanfermann, J. Masso, T. Radke, E. Seidel, and J. Shalf. The 
Cactus framework and toolkit: Design and applications. In Vector and Parallel Processing 
- VECPAR '2002, 5th International Conference, Lecture Notes in Computer Science, Berlin, 
2003. Springer. 

[27] http://www.cactuscode.org. 

[28] Erik Schnetter, Scott H. Hawley, and Ian Hawkc. Evolutions in 3D numerical relativity using 

fixed mesh refinement. Class. Quantum Grav., 21(6):1465-1488, 21 March 2004. 
[29] http://www.carpetcode.org. 

[30] Christian Rcisswig. Binary Black Hole Mergers and Novel Approaches to Gravitational Wave 

Extraction in Numerical Relativity. PhD thesis, Leibniz Universitat Hannover, 2010. 
[31] Denis Pollney, Christian Reisswig, Erik Schnetter, Nils Dorband, and Peter Diener. High 

accuracy binary black hole simulations with an extended wave zone. arXiv.0910.3803, 2009. 
[32] Denis Pollney, Christian Rcisswig, Nils Dorband, Erik Schnetter, and Peter Diener. The 

Asymptotic Falloff of Local Waveform Measurements in Numerical Relativity. Phys. Rev., 

D80:121502, 2009. 

[33] Nigel T. Bishop, Roberto Gomez, Paulo R. Holvorcem, Richard A. Matzner, Philippos 
Papadopoulos, and Jeffrey Winicour. Cauchy-charactcristic matching: A new approach to 
radiation boundary conditions. Phys. Rev. Lett., 76(23) :4303-4306, 3 June 1996. 

[34] Nigel T. Bishop. Linearized solutions of the Einstein equations within a Bondi-Sachs framework, 
and implications for boundary conditions in numerical simulations. Class. Quantum Grav., 
22(12) :2393-2406, 2005. 

[35] B. Szilagyi, Roberto Gomez, N. T. Bishop, and Jeffrey Winicour. Cauchy boundaries in 
linearized gravitational theory. Phys. Rev. D, 62:104006, 2000. 

[36] Miguel Alcubicrrc, Gabrielle Allen, Thomas W. Baumgarte, Carles Bona, David Fiske, Tom 
Goodale, Francisco Siddhartha Guzman, Ian Hawkc, Scott Hawley, Sascha Husa, Michael 
Koppitz, Christiane Lechner, Lee Lindblom, Denis Pollney, David Ridcout, Marcclo Salgado, 
Erik Schnetter, Edward Seidel, Hisa aki Shinkai, Deirdre Shoemaker, Bela Szilagyi, Ryoji 
Takahashi, and Jeffrey Winicour. Towards standard testbeds for numerical relativity. Class. 
Quantum Grav., 21(2):589-613, 2004. 



